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Abstract: This work considers the problem of super-resolution. The goal is to resolve a Dirac distribution from 
knowledge of its discrete, low-pass, Fourier measurements. Classically, such problems have been dealt with parameter 
estimation methods. Recently, it has been shown that convex-optimization based formulations facilitate a continuous 
time solution to the super-resolution problem. Here we treat super-resolution from low-pass measurements in Phase 
Space. The Phase Space transformation parametrically generalizes a number of well known unitary mappings such 
as the Fractional Fourier, Fresnel, Laplace and Fourier transforms. Consequently, our work provides a general super¬ 
resolution strategy which is backward compatible with the usual Fourier domain result. We consider low-pass 
measurements of Dirac distributions in Phase Space and show that the super-resolution problem can be cast as Total 
Variation minimization. Remarkably, even though are setting is quite general, the bounds on the minimum separation 
distance of Dirac distributions is comparable to existing methods. 
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I. Introduction 

Ernst Abbe’s foundational work ||T| in 1873 reported an observation regarding lack of resolvability of optical 
features beyond the diffraction limit. This problem is central to several areas of science and engineering such 
as optics m, 0, imaging a, geophysics 0, depth sensing 0 and astronomy Q. 

In its abstract form, the problem can be stated as follows: How can we recover a iC-sparse signal (spike 
train/Dirac impulses) with unknown locations and amplitudes from the knowledge of its low-pass measurements 
in the Fourier domain? The problem is challenging because it asks for recovery of a non-bandlimited signal from 
its projection onto the subspace of bandlimited functions. Two predominant approaches exist in the literature to 
super-resolution: optimization and parameter estimation based solutions. 
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A: Minimum separation fc'. Cut-off frequency 
Donoho lfT4ll Kahane ifTSll Candes-G. ifTTIl Moitra 1161 





A > ^ 


A > 


1 

fc-l 


TABLE I: Exact Recovery Condition for Super-Resolution. 

■ Optimization Based Super-Resolution: Early roots of £i norm based optimization can be traced back to 
ii. This approach has been revitalized due to the advent of compressed sensing 13. A very recent idea in this 
context is that of continuous sparse modeling where the signal attributes are estimated on a continuously dehned 
grid QOl, O (instead of its discrete counterpart, as is the case for the usual compressive sensing setup). 

Bredies and Pikkarainen ED present an optimization method in the space of signed measures. In their 
setting, they consider discrete measurements while the solution space is inhnite dimensional. In parallel, de 
Castro and Gamboa 03 and Candes and Eernandez-Granda im consider a Total Variation (TV) formulation 
for super-resolution. 

■ Parameter Estimation Based Super-Resolution: The super-resolution problem has been studied in the 
signal processing context with the goal of resolving overlapping echoes/time-delay estimation ini, multi-path 
characterization ||6l and deconvolution. Li and Speed ifTSl considered super-resolution in form of parametric 
deconvolution of spike trains. Vetterli, Blu and co-workers 03, 113 developed the idea of super-resolution 
as a sparse sampling problem. Eldar and co-workers developed super-resolution methods in the context of 
sub-Nyquist sampling and the Xampling framework mi. Einally, since the super-resolution problem is closely 
linked with the spectral estimation problem ifTSl . MUSIC, ESPIRIT and matrix pencil ll22l . l2^ based methods 
have also been used. 

In contrast to parameter estimation techniques where the recovery condition is based on the number of spikes, 
non-parametric methods provide a bound in the form of a minimum separation condition. More precisely, let K 
be the number of spikes to be super-resolved, fc be the cut-off frequency in Eourier domain and let A denote the 
minimum spacing between any two spikes. While parameter estimation methods require fc > 2K + 1 assuming 
that K is known a priori, non-parametric methods super-resolve the spikes provided that fc > /sr (A) where 
/sR is some function. The works of Donoho iflTl . Kahane ITSl . Candes and Eernandez-Granda ITTl and Moitra 
m provide a theoretical guarantee for the super-resolution problem in terms of /sr. We summarize the recovery 
guarantees in terms of /sr in Table |l] 

The Eourier transform is well suited for examining signals which are linear combinations of sinusoids. 
In many practical applications such as radar, sonar, holography, wave-physics and quantum optics, the basic 
building blocks of the signals are not sinusoidal. Often polynomial phase, Eourier-like transformations of form 
gjHt) are well suited for analysis of such signals. Eor this purpose, tools such the Eresnel transform ll24l . 
Eractional Eourier transform ll25l and the Chirp transform 123 . l27l have been developed in the literature. 

Our goal here is to extend existing super-resolution results to integral transformations other than the Eourier 
transform. To this end, we cast the super-resolution problem in phase domain which parametrically generalizes 
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some of the well-known unitary transformations. Some examples are listed in Table [III We show that exact 
recovery of spike trains from low-pass measurements in phase space is possible by minimizing the signal’s 
TV-norm. This is accomplished using a convex program. Remarkably, even with the general construct of the 
problem, our theoretical guarantee for exact recovery is the same as in m. 

Throughout the paper, M, C and Z denote sets of real, complex and integers. We use z* to denote the complex 
conjugate of z G C, ^ (z + z*) is the real part of z and j = '/—I. Discrete sequences are represented by 

s [m] , TO G Z while their continuous counterparts are represented by s (f), f G K. The L 2 inner product between 
functions / and g is denoted by (/, p) = / fg*dt. Function composition is denoted by (/ o g) {x) = fig (x)). 
Convolution between functions / and g is defined as (/ * g) (t) = J f {z) g (t — z) dz. We use s to represent 
the row-vector of a discrete sequence s [to] , to G [0,..., M — Ij. We use bold capitals to represent matrices, 
for example D and is the pseudo-inverse of D. Calligraphic letters such as £ are used to denote operators. 
Sets are denoted by capitalized roman font, S. The estimate of function/sequence /i is represented by Jl. 


II. Phase Space Representation of Signals 


The term Phase Space naturally arises in areas of optics ll28ll and mathematical physics 12911 . The Fractional 
Fourier transform (FrFT) l25l . |[30l-|[32l and the Linear Canonical Transform (LCT) l28ll . |[33], l34ll are 
instantiations of phase space transformations. These transformations have found a number of applications in 
signal processing 1 ^ and communication problems such as shift-invariant sampling/approximation theory IMI . 
l34l . operator theory, multi-carrier communications 1^ . array processing ED and optical signal processing 
l28l . In this paper, we will work with the Linear Canonical Transform (LCT) l2^ . 

Definition 1 (LCT): Let A = [“[]], with ad — bc= 1. The LCT of a function / (f) G K is a parametric, 
unitary, integral mapping, £a ■ f ^ f with respect to the time-frequency kernel k/^, 


f H = -Ca [/] (w) 

LCT 

where the transformation kernel is defined by 


I if, kA 


bfiO, 

5 = 0 , 


kA it,U}) 


1 


V-j27r5 


exp (^~J^ ~ . 


( 1 ) 


( 2 ) 


Now since 5 = 0 amounts to dilating the function f {du) (see ([D), we will develop our results for 

the case when b 0. 

The LCT satisfies a useful operator composition property: £ai o Ca^ = Ca^ with A 3 = A 2 A 1 . This can 
be used to show that. 


fit)=CA-i[f]it) = 

Inverse-LCT 


(^f,kA-i it,-)) 

(ai) 


57^ 0 

5 = 0. 


(3) 


The LCT parametrizes a number of well-known, unitary transformations, some of which are listed in Table HD 
For A = [ _?i 0 ] = Aft, the LCT amounts to the Fourier transform (upto a constant). (FT). Similarly, with the 
2x2 rotation matrix Ag (see Table Hill, we obtain the Fractional Fourier transform (FrFT). Matrix factorization 
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TABLE II: Parametric Representation of Unitary Transformations 


Parameter Matrix (A) 

r COS 9 sin 9 1 _ \ 

L - sin 0 cos 0 J — 

[-to] = ^FT 

[°^]=Alt 

r j cos 9 j sin 9 1 
L J sin 9 —J cos 9 \ 

[l\] 

ivn 

1 r n 


Corresponding Transform 
Fractional Fourier Transform 
Fourier Transform (FT) 
Laplace Transform (LT) 
Fractional Laplace Transform 
Fresnel Transform 
Bilateral Laplace Transform 
Gauss-Weierstrass Transform 
Bargmann Transform 


shows that the LCT is related to the FT and the FrFT. The Fourier matrix factorization is simple: 


a b 
c d 


b 0 
d 6-1 


Aft 


1 0 
a/6 1 


44- A = MiAftM2, 


and leads to the implementation £a = ° ° - The Fractional Fourier matrix factorization is more 

involved. Relying on the Iwasawa Decomposition llT5l . 


a b 
c d 


= 


r 0 

0 r-1 


1 u 
0 1 


44« A = AgDU, 


where D is a diagonal matrix with F = and U is an upper-triangular matrix with u = {ab + cd) /F^. 

A useful operation that is linked with phase space is the convolution/filtering operator. For the Fourier 
domain, we have the convolution-multiplication property: {f * g) — C/^-i [£a [f]^A [ff]], with A = Api- 
Unfortunately, this property is not preserved in phase space. To circumvent this problem, we use a version of 
the FrFT convolution operator 1^ . 

Definition 2 (LCT Convolution/Filtering): Let *a denote the convolution/filtering operation in LCT domain 
and * be the usual Fourier domain convolution operator. Convolution of functions / and g in the LCT domain 
is dehned by 


if*Ag) (t) = (/ (<) * g it) e+^^) . 




(4) 


By following the steps in 
property: 


Convolution of Modulated Functions 

, it is easy to verify that the operation in (|4|i admits a convolution-multiplication 


/^A [/*A5] (w) = e ^^/(w)p(a;). 


(5) 


III. Sparse Signals in Phase Space 
Consider a A-sparse object/spike train modeled by, 

s(t)=y^ Cfe(5(f - 4), f G R (6) 

where 5 denotes the Dirac mass with weights {ck}k G C that is activated on locations {tk} G [0, r), fc = 
0,..., A — 1. Since we are dealing with a hnite length signal s that lives on the interval [0, t), we investigate 
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its representation in phase space using the Fractional Fourier series analog of the LCT. 

Definition 3 (Linear Canonical Series): Let / be a compactly supported function such that / 7 ^ 0, Vf G [0, r) 
and zero elsewhere and let = s/^TTb/r, b 0. The Linear Canonical Series (LCS) expansion of the function 
/ is given by 

n=+oo 2^ 

/ (t) = Kb,T f [n] La {t,nuJob), ujq = —. (7) 

n= —00 

The LCS coefficients / are evaluated at w = nujob, n G Z, 

/ [n] = Kb,rCA [/] {nuJob) = Kb^r{f,kA i^nuJob)). ( 8 ) 


Note that by appropriately parameterizing A, one can easily design the basis functions for any of the 
transformations listed in Table |II] For example, with A = Apr in ®, we obtain the Fourier Series expansion 
off. 

We now compute the series coefficients for the sparse signal s (f) in We use ®, to compute the LCS 
coefficients 


/* 

^[n] = Kb,TCA[f]{nuJob) = Kb,T / s (t) kX{t,nu!ob) 


dt 


.K-l 


E x V — _L 


(9) 


Plugging the coefficients into the LCS representation ([Tl), we obtain the phase space representation of s (t) 

m ^ 27r 

S [t) = 2_^ 5 [mj /CA (i, mccoo), OJQ = — 


m— — oo 


.. at^ m—+oo /K—1 
e J 2 b 






? 7 i=—00 \ A ;=0 


T m— + oo 
e J 2b 


y[m] 


( 10 ) 


( 11 ) 


From lUll . we see that y[m] are precisely the Fourier series coefficients of rs (f) . Thus, even though we 

are dealing with the phase space, the linear frequency modulated, sparse signal s is completely characterized 
by its Fourier series coefficients y. 


IV. Super-Resolution in Phase Space 

A. Problem Formulation 

Let sine (t) = . Consider the frequency modulated function, 

fiiP (t) = (f^/b) sine ((fl/&) t). (12) 

Note that this function is (f27r)-bandlimited because in phase space the function i^lp (t) is compactly supported, 
or, 

ipH = £A[^Lp|H = ^n(^) 










where 11 (w) = 1, |aj| ^1/2 and zero, elsewhere. 
With s defined in (fTTI) . its low-pass version is, 
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h{t) = (s*a0lp) it\ 'S (13) 

Low-Pass Filtering ^ |mKLf2T/26J 

where [-J is the floor operation. 

Suppose we observe N discrete, low-pass measurements sampled with sampling rate T = b/Q,, 

h[n\=h{t)\^^^T^ T = b/n, n = 0,...,N-l. (14) 

By modulating h, we obtain measurements of the form, 

y [n] = \/j2Td)Te~^^ h[n] = y [m] (15) 

Modulated Measurements l"iK/c 

where fc = [nT/ 2 &J. In vector-matrix notation, we have low-pass measurements, y = ViDpiy, where Vidft S 
([ 2 Nx{ 2 fa-i) jg jjjg ^sual inverse-DFT matrix with elements [Vidft]^^ = and we assume that 

N ^ 2/c — 1 so that Vidft is a full-rank marix. 

Having obtained y, the question then is: how many samples of y are sufficient for complete characterization 
of y[m] g-jrnuiotk'} Indeed, from spectral estimation theory fi22\ . we know that we 

need at least 2K + 1 values of y to solve for {ck,tk}k- Consequently, whenever 

(16) 

the system of equations is complete in y meaning that we have at least 2K + 1 values of y and we can solve 
for {ck,tk}- Thus (fTfil l provides a bound on the minimum sampling density in phase space. 


B. Super-Resolution Via Convex Programming 

Given N measurements y, we obtain y using y = V|Qpjy (fTSI) . From the phase space development of the 
problem, we know that 

y [to] '^2ik 0 (17) 

y [t) \m\ ^ fc = li^T/2b\ 

where 

y(f) = ^ S {t - tk) (18) 

Pk 

and Pk '=^ are the new weights for s (t). 

We are now left to solve the standard super-resolution problem in where one has access to the low-pass 
measurements y[m] and the signal to be super-resolved is prescribed in (fTSl l. Hence, the problem of recovering 
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/X from y can now be solved by using. 



Super-Resolution in Phase Space {Primal Problem) 


minllFllTV 

subject to \y[m]^ [ Jl (t) 


1 JO 



(19) 

where, for the model assumed in (fTSI l. the TV-norm amounts to, ||/r||-j-v = Sfe \Pk\ = |cfc|- We note that 
in principle ||/r||jv = ||s||jv, however, due to the nature of phase space measurements, the real/complex weights 
{ck}k need to be demodulated using the linear frequency modulation term which depends on A. 

The problem in (fT^ seeks to recover the infinite-dimensional variable p from finitely many constraints 
set up in dnii. This continuous optimization problem has a tractable dual problem. As was shown in im, a 
semidefinite program (SDP) can be used to recover J1 by computing {tk}k first and then recovering {ck}k using 
a least squares fit. The SDP equivalent im of the convex dual of (fT^ is. 


Semidefinite Program (Duai Problem) 


max 5R (y, u) 

u,M 


M u 

u* 1 


^ 0 , 


subject to. 


y [M] 

,mGSi 


j_ - = Sj 

J 


where Si = [1,2/^ + 1 — j], S 2 = [0,2/^], M S some Hermitian matrix and u S 

is a complex vector. 

The SDP input y results in a vector u. In order to recover the locations {tk}k, we construct the polynomial 
of degree Nq = Afc, 

PNo (z) = 1 - Xl|fe|<2/, ^ ^ 


The roots of {z) ,z = lead to the locations {tk}k- Knowing y together with the estimates, {tk}k, we 
use the constraints in (fTTI) to set up a system of equations which leads to amplitude estimates Ck = 

(see (fT^ l. Finally, we recover our super-resolved signal, 's{t) = ~ ^k)- Stepwise procedure for 

super-resolution in phase space is outlined in Algorithm 1. 

In view of im, let us invoke the definition of minimum distance A = inf{t^}^:t^/ti |ffc ~ ti\- With /c = 
[r2r/2&J , the exact recovery requirement for phase space is as follows: 

Theorem 1 (Exact Recovery in Phase Space): Let the support set of s (t) in ® be S = {tk}k- *^fi® 
minimum distance obeys the bound A (S) fc ^ 2, then s (t) is a unique solution to (fT^ . 

The proof of this theorem is a straight-forward consequence of ca. Moreover, due to inherent Fourier structure 
of the phase space problem, our work may benefit from the ideas discussed in o, na, di, ms. 


C. Remarks and Discussion 

■ Backward Compatibility With the choice of parameter matrix A = Apr (cf. Table in|), our result coincides 
with the usual, Fourier domain case of super-resolution ini, m. Furthermore, A = Ag relates to the case of 
Fractional Fourier domain for which our result generalizes a previous known result llJTlI . 
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Algorithm 1: Super-Resolution in Phase Space. 


Input: Low-pass samples h [n] = {s 0lp) 

_ , a(nT)^ 

Modulate Samples h [n] —)• yjj2'Kb |r|e''''^ h[n] = y [n] 

Data: y = y € 

Solve SDP : max 5ft (y,u) subject to ^ V ^ 

u,M [ U i 

Construct Polynomial: piVo ^ ~ ^\k\<2fc 

Support: pNg = q ->■ 

Weights: nim\y[m]-J 2 k ^0 I 

Pk I I 

Output: s{t) = Cfc(5 [t-tk), {ck = 


■ Exact Recovery Condition Even though our super-resolution naturally extends to a number of well known 
unitary transformations, the exact recovery condition remains unchanged. Hence re-formulating the super¬ 
resolution problem in context of phase space comes at no extra cost in the sense of recovery requirement. 


D. An Application of Super-Resolution in Phase Space 


Bandlimted signals are compactly supported in the Fourier domain. When a bandlimited signal is corrupted 
by additive impulsive noise or AIN, the holes/zeros in the spectrum are filled by the spectral components that 
characterize the impulsive noise which is essentially non-bandlimited. Wolf IMl used the idea of curve-fitting 
the out-of-band components for identification of impulsive noise components. Here, we formulate the problem 
of denoising linear frequency modulated (LEM) signals that are corrupted by AIN. Since LEM signals are the 
basis functions of phase space transformations, it is clear that such signals are bandlimited in the LCT domain. 

Consider a bandlimited LEM signal 


rsL (t) = Kb,T y^, , ^^BL M kA {t, muiQb), 

with Ebl [™] = 0, \m\ > M and let r (t) = pbl (t) + s (t) be the signal corrupted by AIN. Clearly, r (f) is 
non-bandlimited in phase space due to s (t). Suppose we observe low-pass filtered samples of r (f), that is. 


_ a{nT)^ 

{r*A4'ip) {nT) = e ^ 

V ^ > 

h[n\ 


(Ciyi [to] -f 02^2 M) 

I ^ /c r 1 

Vr[m\ 


where ci,C 2 are known constants, yi [to] = tbl [to] e~^' 
[Hr/26J = [r/2Tj. Again, let us define y[n\ = h [n] 
fc>M + 2K + 1, we have. 


2 ° and y 2 [to] = y[m] (as in (fTTl) ) where fc = 
^ 2 b' ,n = 0, l,iV> 2/c + 1. Provided that 


{ ciyi [to] -f C2y2 [to] [to] < M 
02% [to] [to] > M 

which leads to complete characterization of s (t) since the 2K -f 1 values of % [to] can be used with ( [19] ) to 
solve for s (t). With y = y^, to > M we can use Algorithm 1 for exact denoising of r (t). 
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V. Conclusion 

We develop a method for super-resolution in phase space. The phase space transformation generalizes a 
number of well known transforms (see Table HU). More precisely, we are concerned with recovery of spike trains 
from their low-pass samples. For this purpose, we filter the spike train with a kernel which is bandlimited in phase 
space. We show that even though we are dealing with a general class of parametric transformations, the low- 
pass samples are completely characterized by chirp-modulated Fourier series. Having made this link, we show 
that the recovery of spikes from their low-rate measurements can be cast as a total-variation minimization—a 
problem that can be tackled by convex programming. In closing, our work extends the recent results of ca 
without altering the exact recovery condition. That said, the cut-off frequency is a function of the transform 
being used for investigation. Our work warrants future research, specially for the case of additive noise. 
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